% Select TF dataset (containing indices indicating the row numbers of TF
% target genes in the expression data file
TF = Znclustarg;

% Prepare data for selected TF
% Select corresponding expression data from expression file containing TFs
% in the rows and conditions in the columns
TFdata = allexprdata(TF, :);
TFdata = TFdata';

% Find correlation between expression patterns of all pairs of TFs
AllTFcorr = corr(TFdata, 'rows', 'complete');
figure; h = pcolor(AllTFcorr);
set(h, 'EdgeAlpha', 0);

% combine all correlations into one vector
concTF = [];
for i = 1:length(AllTFcorr(:,1))
concTF = [concTF AllTFcorr(i, :)];
end
%remove self correlations from this matrix and sort
concTF = concTF(find(concTF < 1));
concTF = sort(concTF);
figure; hist(concTF,30);

% Calculate correlation at top 5% of all TFs and enter result here
D = 0.3069;

% Find the number of pairs in test set with correlation greater than this cutoff
p = length(find(concTF > 0.3069))

% Expression coherence is the pairs above threshold / total pairs
EC = p/length(concTF)

%%%% For significance: testing of random sets of same size
% Test 1000 random sets and calculate EC for each for current set
% size
setsizes = [4; 5; 7; 8; 9; 10; 12; 14; 17; 27; 38; 68; 133];

EC = zeros(1000, length(setsizes));

count = 0;
for s = 1:length(setsizes)
    size = setsizes(s)

    count = count+1;
    for trial = 1:1000
        trial
        TF = round(rand(size, 1).*(length(allexprdata(:,1))-1))+1;
        % Prepare data for selected TF
        % Select corresponding expression data
        TFdata = allexprdata(TF, :);
        TFdata = TFdata';

        % Find correlation between expression patterns of all pairs of TFs
        AllTFcorr = corr(TFdata, 'rows', 'complete');
     
        % combine all correlations into one vector
        concTF = [];
        for i = 1:length(AllTFcorr(:,1))
            concTF = [concTF AllTFcorr(i, :)];
        end
        %remove self correlations from this matrix and sort
        concTF = concTF(find(concTF < 1));
        concTF = sort(concTF);
   
        % Calculate correlation at top 5% of all TFs and input here:
        D = 0.3069;

        % Find the number of pairs in test set with correlation greater
        % than this cutoff
        p = length(find(concTF > 0.3069));

        % Expression coherence is the pairs above threshold / total pairs
        % EC(trial, count) = p/length(concTF);
        EC(trial) = p/length(concTF);
    end
end